Nonlinear electronic transport in nanoscopic devices: 
Nonequilibrium Green's functions versus scattering approach 
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We study the nonlinear elastic quantum electronic transport properties of nanoscopic devices 
using the Nonequilibrium Green's function (NEGF) method. The Green's function method allows 
us to expand the I — V characteristics of a given device to arbitrary powers of the applied volt- 
ages. By doing so, we are able to relate the NEGF method to the scattering approach, showing 
their similarities and differences and calculate the conductance coefficients to arbitrary order. We 
demonstrate that the electronic current given by NEGF is gauge invariant to all orders in powers 
of V , and discuss the requirements for gauge invariance in the standard Density Functional Theory 
(DFT) implementations in molecular electronics. We also analyze the symmetries of the nonlinear 
conductance coefficients with respect to a magnetic field inversion and the violation of the Onsager 
reciprocity relations with increasing source-drain bias. 

PACS numbers: 72.10.-d,73.23.-b,73.63.-b,85.65.+h 



I. INTRODUCTION 

There is a growing experimental and theoretical in- 
terest in quantum nonlinear electronic transport prop- 
erties of nanoscopic devices. Experiments in mesoscopic 
semiconductors, such as quantum dotsjii^ and quantum 
rings^£ have focused on the investigation of rectifica- 
tion effects and violations of the Onsager-Casimir reci- 
procity relations. Similar issues have also been examined 
in the electronic transport through carbon nanotubes. 6 
Nonlinear transport is also of major interest in molecular 
electronics j-£ where considerable experimental effort has 
been put in showing that single molecules can be used as 
diodes, transistors, and switches. 

The development of a comprehensive quantum nonlin- 
ear electronic transport theory, a non-equilibrium quan- 
tum many-electron problem, is still quite a challenge. 
Notwithstandingly, significant advances have been al- 
ready achieved, particularly by restricting the theoret- 
ical analysis to elastic processes. Within this approx- 
imation, theoretical progress has mainly been achieved 
by pursuing two apparently very different paths, namely, 
the scattering approach put forward by Biittiker and 
collaborator s 8 ' 9 ' 10 : 11 ' 12 and the nonequilibrium Green's 
function methodJ^i^±£ 

In both approaches the current-voltage I — V charac- 
teristics is written in terms of transmission coefficients 
that account for the potential landscape built in the de- 
vices due to the applied bias. The scattering approach 8,9 
casts the current as a power series of the bias. The S- 
matrix serves not only to compute the transmission, like 
in the Landauer formula, but also to calculate the elec- 
trostatic potential built in the conductor by means of 
physical considerations. Those guarantee that the elec- 
trostatic potential is gauge invariant order by order in 
powers of V. Alternatively, NEGF has also been used to 
investigate the linear and non-linear transport properties 



of mesoscopi o 13 : 15 ' 16 and molecular systemsi 14 ' 17 ' 18 ! 19 : 20 
Here, the current and the electrostatic potential are cal- 
culated self-consistently to all orders at once. The for- 
malism is quite powerful and robust but, as it is often the 
case in self-consistent calculations, different physical pro- 
cesses become inextricable making difficult to understand 
their role and importance for the electronic transport. 

These considerations raise a natural question: To what 
extent are these approaches similar? One of the main 
purposes of this paper is to answer this question and 
explicitly show that both approaches are, in principle, 
equivalent. Furthermore, we show that differences ap- 
pear depending on how the underlying many-body elec- 
tronic problem is approximated. 

To this end we use NEGF to write the current flowing 
through a multi-lead elastic conductor as a power series 
of the bias V, as done in Ref. [l3|. Treating the many- 
electron problem in the Hartree approximation we ex- 
plicitly show that, in the Thomas- Fermi limit, the leading 
nonlinear correction in the I — V characteristics reduces, 
almost exactly, to the scattering approach result. We 
discuss gauge invariance and the Onsager-Casimir reci- 
procity relations. At every step, we also analyze these 
symmetries beyond the Hartree term by addressing the 
standard Density Functional Theory (DFT) implementa- 
tion for molecular electronics^ We stress that, in con- 
trast to the scattering approach, the NEGF formalism is 
not restricted to a local approximation. The many-body 
nature of NEGF allows one to extend its standard imple- 
mentation in a variety of ways. To show this, we address 
the case where the many-body problem is treated in the 
Hartree-Fock approximation.— 

The presentation of the paper is organized as follows: 
In Sec. |H] we present the scattering approach highlight- 
ing its main elements and results for later comparison 
with the NEGF approach. In Sec. IIIII we present the 
model Hamiltonian considered in this study. In Sec. IIVI 
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we calculate the transmission coefficients and the electro- 
static potential using noncquilibrium Green's functions. 
By these means, it is possible to systematically calculate 
the conductance coefficients and the characteristic po- 
tentials self-consistently to all orders in V, as discussed 
in Sec. [V] The similarities between both approaches are 
discussed in Sec. IVI1 where we also show how differences 
appear. The NEGF implementation of the Hartree-Fock 
approximation is described in Sec. IVIIl Finally, our con- 
clusions are presented in Sec. IVfffl 



II. THE SCATTERING APPROACH 

Let us consider a conductor connected by leads to 
a = 1, • • • , N electronic reservoirs at a temperature T. In 
the absence of an applied bias, the system is in thermal 
and chemical equilibrium, characterized by a chemical 
potential fj,g. By applying voltages {V a } to the reservoirs, 
the system is driven out of equilibrium and an electronic 
current flows. 

According to Biittikcr, - in the absence of inelastic pro- 
cesses, the current at the lead a is given by 



where u a p...(r) is the characteristic potential defined by 
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dEfp{E)A a p{E,U{v)) . (1) 



Here fp{E) = fo(E-eVp), where f (E) = (e^^+l)- 1 
is the Fermi distribution function and fee is the Boltz- 
mann constant. For notational convenience, we consider 
Vp as measured with respect to the equilibrium potential 
fx , namely, Vp —> Vp - (M>/e. 

The transmission is identified with 



A a p(E,U(r)) = Tv[l a S a p 



(2) 



where S a p(E,U(r)) denotes the scattering matrix with 
lines and rows associated with the transversal modes at 
the contact a and (3, respectively. l a is the identity 
matrix whose rank is given by the number of propagating 
channels in the contact a. The trace runs over all open 
channels in a and (3. 

The transmission coefficient A a p and the scattering 
matrix S a p are functions of the electron energy and func- 
tionals of the electrostatic potential U (r) in the conduc- 
tor. In linear response, A a p is computed at the equilib- 
rium potential U eq {r) that is established when all reser- 
voirs have the same chemical potential fiQ. Beyond this 
regime, it is necessary to compute U(r) self-consistently, 
as pointed out by Landauer^ 

To make analytical progress, it is convenient to expand 
all quantities in powers of V. The local electrostatic po- 
tential U(r) reads 

U(r) = U eq (r) +J2u a (*)V a 

Qi 

+ ^J2 U M V ^ + 0(V 3 ) (3) 

a/3 



i a p...(r) 



d d 



dV a dv. 







U(r) 



{v y }=o 



(4) 



Here {V^} = is a shorthand for Vy = 0, for all 7. 

Some properties of the characteristic potentials fol- 
low directly from simple physical considerations. For in- 
stance, u a (r) has the following properties^ (a) Changes 
in the electro-chemical potential of the reservoir a should 
not affect U (r) inside f3, hence, u a (r) = 0, when r is taken 
inside the reservoir f3 ^ a. (b) For r inside the reser- 
voir a, U(r) = V a and, thus, u a (r) — 1. (c) A global 
change of the applied potentials, V a — ► V a + Vq, makes 
U(r) — » U(r) +V~o, implying the sum rule J2 a u a( r ) = 1 
for all r. 

The current I a , written as a power series of the applied 
voltages, is cast as a function of the coefficients Q a p—, 
namely, 



I a = / ]GapVp + /lG a p-yVpVy- 
(3 n 



■E 

0jS 



GapysVpVyVs- 



(5) 

In line with the standard notation^ we do not write I a as 
a Taylor series in {V^J. In Sec. \V\ we will see how such 
notation determines the symmetrization of the indices 
ot,/3, • ■ ■ in the conductance coefficients Q a p.... 

The coefficient Q a p corresponds to the linear conduc- 
tance, as given by the Landauer formula 



Gu0 



h 



dE 



A aP (E,U cq (r)). (6) 



Here A a p(E, U cq (r)) is the multi-lead Landauer-Buttiker 
transmission coefficient, 8 with S computed using U cq (r), 
as standard. 

The first non-linear current correction, represented by 
Ga0i reads&2 



Qa0i — —r~ 



2e 3 
h 



dr 



u 7( r ) " 



SA a , 



eSU(r) 



{K,}=o 



, (7) 



where the spatial integration is taken over the region 
where 6A a p/5U(r)\ry a y =0 is non vanishing, namely, in- 
side the conductor. 

The above expression depends explicitly on the elec- 
trostatic potential via u a (r). To determine u a (r), the 
formalism has to be supplemented by a self-consistent 
microscopic electronic structure calculation, or by an ad- 
equate approximation. The latter was constructed in 
Ref. H by using the following argument: The potential 
U(r) is related to the bias generated electronic density 
imbalance Sn(r) in the conductor. In turn, Sn(r) arises 
from the charge injected by the leads and the induced 
charge in the conductor, in response to the injected one. 
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The injection properties of the sample are given by the 
injectivity, which reads 



drf(r, a) 
dE 



1 
2~7ri 



x£Tr 





dE 



dfo 
dE 



'/3a 



eSU(r) eSU(r) 



SS/3 a 



(8) 



evaluated at {V^.} = 0. The superscript s labels quan- 
tities obtained within the scattering approach. Since 
dn s /dE includes spin degeneracy, our definition differs 
from the one of Ref. d by a factor 2. 

To linear order in V , the induced charge density is 
given by 



dn ind (r) 



J dr / n(r,r')u a (r')dV a 



(9) 



where IT(r,r') is the Lindhard polarization function. 23 
The scattering approach does not provide a recipe to ob- 
tain the later. However, by recalling the relation between 
the Wigncr-Smith time delay and the conductor density 
of states, c£nj n d(r) can be readily written in the Thomas- 
Fermi approximation as 



dn ind (r) = e 



dn s (r) 
dE 



u a (v)dV a . 



The local density of states dn s /dE is 



dn s (r) 
dE 



E 



rfn s (/3,r) 
dE 



(10) 



(11) 



where dn s ((3,r)/dE is called emissivity and is given by 



dn s {(3,v) _ _L_ 



dE 



dE - 



5> 



dE 



P a e8U(r) 



Ml* 
eSU(r] 



These elements render the Poisson equation 



V 2 w Q (r) 



47re 



2 dn s (r, a) 
~~dE ' 



(12) 



(13) 



where both the density of states and the injectivity de- 
pend only on the scattering matrix. 

Higher order conductance coefficients Q also 
be calculated in a straightforward way. Obtaining 
self-consistent equations for the characteristic potentials 
u a f3~y... becomes increasingly more involved, but is still 
possible within the Thomas-Fermi approximation. 24 

In the next Sections we define a general model for a 
conductor and use the NEGF approach to show how 
to systematically obtain the coefficients Q a gj... and the 
characteristic potentials u Q| g 7 ... to arbitrary order. 



III. MODEL HAMILTONIAN 

We separate the system in two regions, namely, the 
leads (L) and the conductor (C) to write the Hamiltonian 
as 



TL = Tih + Tic + Wlc ■ 



(14) 



For dcfinitcncss, we introduce a surface S enclosing the 
conductor, to partition the model Hubert space. The 
lead Hamiltonian reads 



s kaas kaas ' 



(15) 



where k is the electron transversal wave number at the 
channel a (a — 1, • • • , N a ) in the lead a (a — 1, ■ ■ ■ , N). 
The electron spin is s =T, J. and the c\ aas {c kaas ) are the 
usual fcrmionic creation (annihilation) operators, with 

{ c Las' c k> a >a>s>} = 6kk>6 aa 'S aa ,5 ss/ . The threshold en- 
ergy to open the transversal propagation mode a in the 
lead a is E aas . We assume free motion in the direc- 
tion along the leads. Hence, E^aas = E aas + h 2 k 2 /2m*, 
where m* is the electron effective mass. The electrons at 
the lead a are in thermal equilibrium with the reservoir 
at temperature T to which the lead is connected. This 
reservoir is characterized by a chemical potential /i Q . 
The conductor Hamiltonian reads 



(16) 



where dj ls (d IJlS ) creates (annihilates) an electron at the 
/x-th state of an arbitrary basis {v} that spans the con- 
ductor eigenstates. Since we consider He as a bilinear op- 
erator, electron-electron interactions are only taken into 
account in the mean-field level. This is a good approxi- 
mation, provided the system is open^ and hence neither 
charging nor electronic correlations effects are expected 
to play an important role. 

The term that couples the leads to the conductor is 



H-LC 



kaas^fJ.s 



H.c. 



(17) 



When there is a difference between the reservoirs' 
electro-chemical potentials the system is driven out of 
equilibrium and a current flows. In the stationary regime 
a time-independent non-equilibrium self-consistent elec- 
trostatic potential U(r) is formed. It depends on the 
applied bias, as well as on the system geometry and ma- 
terial properties. In most of the paper, we assume U(r) 
to be local. A non-local U is discussed in Section IVlIl 

The Hamiltonian Wl of Eq. (fTS"]) assumes free propa- 
gation in the leads. Alternatively, without significant in- 
crease in complexity, it can also represent periodic semi- 
infinite leads In any of these events, the surface sep- 
arating conductor and leads has to be chosen in such a 
way that, at any given lead a, U(r) « V a . As a con- 
sequence, the spatial dependence of U(r) is entirely ac- 
counted for by Tic- This construction not only limits 
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the arbitrariness in defining the model Hilbert space, but 
also guarantees a simple prescription for computing the 
characteristic potentials u a 0...(r), as we shall discuss in 
Section [V] In the standard DFT approach for molecu- 
lar electronics, although not emphasized, this is the key 
notion behind defining an "extended molecule" and it is 
essential to ensure gauge invariance.— 

There is a more basic principle behind the above 
construction® than just simplifying calculations: The sur- 
face S defines the leads region at a position where no elec- 
trical field lines penetrate its surface. Hence, the charge 
within the volume V enclosed by iS is constant. 

For the sake of simplicity we restrict our considerations 
to weak magnetic fields, or more precisely, to systems 
where we can neglect spin-orbit and Zeeman interactions 
and consider only orbital effects due to an external mag- 
netic field. Hence, in what follows, except for Sec. IVII1 
we omit the spin index and replace the sums over spin 
projections by their degeneracy factors. 

We also do not explicitly include the possibility of 
capacitive couplings in our model. In pumping experi- 
ments, such kind of coupling is likely to dominate the 
transport,— as discussed in Refs. [28|. In dc nonlinear 
transport, the effect of setting a fixed back gate voltage 
and, hence, defining bias mode, can be relevant in non- 
linear conductance of quantum dotsJ^ We stress that 
both the scattering and the NEGF approaches can eas- 
ily accommodate situations where the number of leads 
does not coincide with the number of gate voltages {V^} 
and/or devices with more than one conductor. 



IV. THE NEGF APPROACH 

For elastic processes, the electronic current at the 
leads can be written in terms of the conductor Green's 
functions^ namely 



r(a) 



la = -'- 



2e 



— ImTr(r Q [G<(E) + f a (E)G r (E)] }, 

(18) 

where symbols in bold face correspond to matrices whose 
rows and columns are states of the conductor basis set 
{fi}. The Green's functions G r (E) and G<{E) are 
the Fourier transforms of G r „ v (t — t') = —{i/Ti)6(t — 

0<{d„(*),4(* , )}> and G< v (t-f) = (i/n)<4(fX(t)) 
respectively, where (• • ■) is defined as standard^ The 
decay width or line width matrix elements are given by 



In the energy representation, the Green's function G^ 
is given by 



G r{a \E) = \eI - H c - £ r(a) (£) 



(20) 



The conductor lesser Green's function G<„ follows di- 
rectly from the Dyson equation^ 



G<{E) = G r (E)'E < (E)G a (E) 



(21) 



For the sake of definitcness the basis set is truncated 
and the matrices have rank M. I is the identity matrix. 
(Later on we shall also use conductor Green's functions 
in the coordinate representation and drop the boldface 
notation.) 

The self-energy matrix elements read 



^fj,,A E ) = E Vfrkaa 9kaa( E ) ^v,kaa 



(22) 



kaa 



where E r ( a ) (and E < ) are obtained by identifying the free 
electron Green's function in the leads g kaa , with 



9La(E) = if a (E)S(E - E kaa ) 



and 



9 r k { :l(E)=Tin6(E-E kaa )+PV 



1 



E — E kaa ' 



(23) 
(24) 



where PV stands for principal value integral. 

The coupling matrix elements V^.kaa are, in general, 
smooth functions of the wave number k and, hence, of 
Ekaa- Using (|24[) and assuming that p aa (E) has a broad 
band width and a smooth energy dependence, the matrix 
elements become energy independent and read 



(25) 



For situations where the broad and flat band approxima- 
tion does not hold, the results we obtain for the I — V 
characteristics have to be modified in a straightforward 
way, as indicated later on. 

In analogy, the self-energy matrix elements E<„ are 
given by 



^u{E) - ^2^2V^Maa9kaa( E )K,kaa 
a — 1 ka 
N 

^EMp- (28) 



I r «U = 27r E V »M« Pk aa {E) V* kaa , (19) 

where p aa (E) is the density of states of mode a at the a 
contact. 

Since the conductor Hamiltonian TCq is a bilinear op- 
erator, the exact conductor Green's functions can be ob- 
tained in closed form using, for instance, the equations- 
of- motion method; ?9 (For a recent review, see Ref. l26h 



Within the wide-band approximation, one arrives to 

V<(E)*if a (E)T a . (27) 

Inserting G r from Eq. (f2"0"|) and G< from Eq. (|2"Tj) into 
Eq. |[T8)) . we write the current I a as in Eq. |T]), 



I a = 



9 ' r°° 

-E / dE fp{E) T a p(E, {V 7 }) 



(28) 
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(29) 



with transmission coefficients given by 

T a/3 (£,{^ 7 }) = Tv[r a G r (£0(r<W-i>)G°(£0 

By means of the useful relation 

G a G r = iGTG' = iG r TG a , (30) 

we obtain X^=i T a (3 = 0, and show that Eq. (|2"8|) satisfies 
current conservation, Yla=tl& = 0- 

As pointed out in Ref. [Tj, the current I a given by 
Eq. (|28[) is invariant under a global shift of the potential, 
that is, V-y -> F 7 + Vo for all 7's and [7 -> U + V . 
This is not sufficient to prove that the formalism is gauge 
invariant. We still have to show that the same condition 
holds for the electron density n(r). This is done in what 
follows. 

The applied voltages {V^.} control the conductor 
charge distribution n(r) and the electrostatic potential 
U(r). The latter, in turn, enters the calculation of the 
conductor Green's functions. Both quantities, n(r) and 
G < , are related by 

n(r,t) = ]T (V4(r,t)^(r,t)) = -2<ft(r|G<(i,t)|r) , 

( 31 ) 

where the factor 2 account for the spin degeneracy. 
By taking ip s (r, t) inside the conductor as ip s (r, t) — 
d^ s (t)(fj,\r), we obtain n(r,t) in matrix representa- 
tion, namely, 



n(r,i) = -2^(r| M >G<„(M)Hr>. 



(32) 



For stationary processes, where G < (i,t) = G < (0), the 
electronic density becomes 



n(r) 



/°° rlE 
^—(r\G<(E)\r) 



(33) 



with obvious matrix representation. In order to close 
the calculational procedure, a relation between n and U 
is needed. This can be done at different approximation 
levels. 

In the Hartree approximation, the electronic density 
n(r) and U(r) are related by 

V 2 U(r) = -4ire n(r) 

(r\G<{E)\v), (34) 

-00 27T 

The boundary conditions are obtained by recalling that, 
by construction, U(r) is constant outside the conductor: 
U(r) = V a when r is taken at the lead a. In addition, the 
problem must to be solved self-consistently, namely, U (r) 
enters the Hamiltonian Tic, which determines G < (E), 
that in turn gives U(r). 

Equation (|3"3"|) leads to a local description of the elec- 
trostatic potential. In Section fVIIi we show how NEGF 



deals with a non-local potential due to the exchange in- 
teraction. 

Equation (|33|) also plays a key role in the standard 
implementations of the density functional theory (DFT) 
in molecular electronics (see, for instance, Ref. [lj for 
a review). In DFT, U(r) = U[n(r)] is considered as a 
functional of n(r), containing exchange and correlation 
interactions in addition to the Hartree one. Accordingly, 
the single-particle states {/i} become Kohn-Sham orbital 
states. Although this is a very appealing construction, 
it is not as sound, from the conceptual point of view, 
as the derivation presented here: DFT is not a mean- 
field theory and a bilinear Hamiltonian, like Eq. (fT4|) . is 
not one of its underpinning elements. A good discussion 
about the shortcomings of the standard DFT approach 
to conductance can be found in Ref. [3]]. 

We conclude this Section by stressing that in local ap- 
proximation schemes Eq. (|33[) is manifestly gauge invari- 
ant, provided the partition given by Eq. (fT^|) satisfies 
the conditions discussed in Section HTT1 In this case, any 
global voltage shift can be absorbed by the energy in- 
tegration, provided U(r) is calculated self-consistently. 
For systems where electronic correlations are build across 
the partition S, gauge invariance calls for a more careful 
analysis. This is the case, for instance, in Kondo systems. 
In the mean field limit, discussed here, such correlations 
are absent. In DFT-NEGF, any semi-local functional of 
the exchange and correlation functional, U xc [n(r)], al- 
lows for a partition S and hence pre serves gauge invari- 
ance, as nicely discussed in Ref. |3li Nonlocal interac- 
tions present in the exact XC functional jeopardize the 
partition construction and spoil gauge invariance due to 
an XC contribution to the characteristic potentials in 
the contacts! 31 ! 33 In such situations, apartition free ap- 
proach, like the one discussed in Ref. |32j is more suited. 



V. LINEAR AND NONLINEAR 
CONDUCTANCE COEFFICIENTS 

In this Section we present a systematic approach to 
calculate the conductance coefficients Q a p-y— and discuss 
some of their properties. For that purpose, we expand 
both f a (E) and T a(3 (E,U(r)), in Eq. (^SJ), as powers of 
the voltages {V a }. 

We start writing the retarded (advanced) Green's func- 
tion as 

G r{a) {E) = — - (35) 

without choosing a particular representation. Next we 
expand G r ^ in terms of the differences between the non- 
equilibrium and equilibrium U and S. As a result, we 
obtain the Dyson equation 



G r(a) = G r(a) Q^y Qr(a) 



(36) 



where the effective perturbation potential V e s is given by 



V cS (E) 



ell — ell ei 



6 



£ £ [x#»>(25 - eV a ) - K {a) (E)] , (37) 



and the equilibrium Green's function by 



1 



E-H Q -eU cq -Y? cq a) {E) 



We now proceed by writing the self-energy 



(38) 



(39) 



as 



E^)(£-eV Q ) = K (a) (E)-eV a 



9Sa 



r(o) 



+e z K 



2t/2 a 



dE 



v a =o 



v a =o 



(40) 



where ^(^(E) = T, r cq a) (E). From Eq. (gDJ) we see that, 
in the flat band approximation, retarded and advanced 
self-energies do not depend on {V a }, since 



r(a) 



dE 



v a =o 



idT a 
2 dE 



. 



(41) 



Hence, V e g depends only on the characteristic potentials, 
namely 



V e g = e^2 u a V a + -e ^2 UapVaVp 



(42) 



a/3 



Inserting the above expression into (|36| we formally ob- 
tain G r( - a ) to arbitrary order in V. 

We now turn our attention to the electronic density 
imbalance 6n(r) = n(r) — n eq (r), that ultimately allows 
us to calculate the characteristic potentials u a g.... To 
obtain Sn(r), we expand G < = G r S < G a in powers of 
{V a } taking, as above, the wide flat band limit. In this 
approximation, the lesser self-energy S < , Eq. (|2o| . reads 

= iJ2fo(E-eV a )T a 

a 

Finally, G < reads 
G< = if G r TG a 



(43) 



-e£V a 
-0{V 2 ) . 



.5/o, 



l ~0^G r T + /o (GqUq,Gq — GqUq,Gq) 



(44) 



Close to equilibrium, when {V a } — > 0, G < 
— 2i/oImGQ, as given by the fluctuation-dissipation 



theorem^ We use Eqs. (|33|) and (|44]) to write the elec- 
tronic density as 



»W = ^nW(r) 



(45) 



where the ts stand for the implicit powers of V 1 . Note 
that n(°)(r) = n cq (r). 

We are now ready to identify the conductance coeffi- 
cients Q a p-f--- order by order: By plugging Eqs. (|36[) and 
(|42p into (|29p we obtain the transmission T a p in terms of 
equilibrium Green's functions and the characteristic po- 
tentials. The later can be computed from G < , as given 
by Eq. (|4~4"|) , with the help, for instance, of the Hartree 
equation 



A. Linear conductance coefficients 

To linear order in V, the current at the contact a is 

(46) 



4 1) = £ft^ 





The conductance coefficients Q a p are obtained from the 
linear expansion of the current (f2"8"]l in the applied volt- 
ages V a . They read 

Q aP = ~ jT d£ T a/3 (£, {F 7 } = 0) . (47) 



where 

T a /3(E, {Vj} = 0) = Tr[r a GS(£)(r<5 a/3 - r /3 )Gg( J B) 

which is identical to the multi-lead Landauer-Buttikcr 
formula^ as it can be verified following, for instance, 
the path presented in Ref. 

Owing to physical considerations the linear conduc- 
tance coefficients Q a p follow some simple sum rules£ Cur- 
rent conservation implies that ^ Q Q a p = 0. This sum 
rule is automatically satisfied by Eq. (T4"5)) , since Eq. (j2"5|) 
does it to all orders in V. Current invariance under a 
global voltage shift V a — > V a + Vb leads to Q a p = 0. 
This is fulfilled, since Y,p T a p{E, {Vy}) = 0. 

Equations (f3"3")) and (j4"4")l give the electronic density as 



dEf (E)lm(r\G r (E)\r) 



(49) 



It is worth remarking that, even in the Hartree approx- 
imation, depending on the system, the computation of 
U cq (r) can already be a formidable computational task. 
In mesoscopic physics, due to the chaotic and/or weakly 
disordered nature of the addressed systems, quantitative 
results can be obtained by a statistical treatment us- 
ing random matrix theory or diagrammatic techniques. 
In molecular electronics a full electronic structure calcu- 
lation is already necessary. 
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An important symmetry of linear transport is unveiled 
by considering an external magnetic field. The linear 
conductance coefficients fulfill the Onsager-Casimir reci- 
procity relations under magnetic field inversion* 3 ^ This is 
indeed the case of Eq. ([48]) i 34 i 37 Using "microreversibil- 
ity" 

G r {a \-B) = [G^\ B )f (50) 

and the cyclic properties of the trace in l|48p. one can 
show that 

G a p{~B) = g Pa (B) (51) 

for a / p. For the diagonal coefficients, where a = /3, one 
can either use current conservation and (|5T|) . or directly 
use Eq. ([30]) to show that G aa (B) = G<xa{—B). 

In the two-terminal case, the conductance itself is an 
even function of the applied magnetic field, namely, 

Qx2{-B)=g X2 {B), (52) 

which is a more stringent symmetry than the reciprocity 
relation for the general multi-lead case. In the scattering 



approach, Eq. ([52]) can be viewed as a consequence of the 
S'-matrix unitarity , 34 ' 37 whereas using NEGF it follows 
from current conservation and from Q aa {B) = Q a a{—B). 
The even symmetry of Q with respect to B has been ex- 
perimentally established 38-39 and has important implica- 
tions for interference experiments in two-terminal meso- 
scopic rings: It does not allow one to measure phase dif- 
ferences, since in the observed Aharonov-Bohm conduc- 
tance oscillations the phase shift is locked either to or 
to n^JSL 



B. Second order terms 

The next order current term, in powers of V, is 

l£> =Y, Q °f* V f> V T (53) 

The coefficient Gap-y, obtained by expanding ^2pfpT a p 
of Eq. ([28]). is formally given by Eq. 0: 



Ga/3j 



2e 3 



dE 



dfo 
dE 



dr 



v 



u 7 (r) 



ST, 



a 6 



e5U(r) 



{V*}=0 



Its explicit expression in terms of equilibrium Green's functions is^ 3 - 
2e 3 r ,„( dfo 



G 



Pi 



dE 



dE 



Tr T a Gl 



G r (TS aP - I» + (T6 afJ - I»Gg 



r 



Ga 




(54) 



(55) 



Note that Ref. [l3| accounts for an energy dependence in 
the self-energy and presents a slightly more general equa- 
tion for Gaf3"f- 

Current conservation implies that ^2 a G a p-y = 0. It 
is straightforward to verify that (p)5"|) satisfies this sum 
rule. The current invariance under a global shift of 
the applied voltages gives a second sum rule,— namely, 
^ 7 (C/ Q( 3 7 + Ga-yp) = 0. We show that the coefficients 
G a f3-f of Eq. ([55]) also respect this sum rule, provided 
that ^2 ti 7 (r) = 1. The latter is also be directly inferred 
from gauge invariance, see Section ITT1 

We now turn our attention to the electron density 
n^(r) of ([45]) . The expansion of G < in powers of V 
and ([55]) give 



1 Hv) = eY j V a \f rfr'n(r,r'K(r') 
Jv 



dn(r, a) 
dE 



(56) 



where 

n(r,r') = -2i 



dE 



fo 



[r\G r (E)\r')(r'\G r Q (E)\r)-R.c. 



(57) 



is formally identified with the Lindhard function and 
dn(r, a) 



dE 



° ^ f-H Hr|G5r a Gg|r) ,r„s) 



2?r 



dE 



is a partial local density of states, called injectivity in the 
scattering approach. The compact form of ([B"5|) is due to 
the flat and wide band approximation. As discussed in 
Ref. [HI, dn(r,a)/dE can be easily modified to account 
for a system specific energy dependence of self-energy S. 
Such corrections are potentially important in molecular 
electronics, where the details of the contacts should mat- 
ter. 

In the Hartree approximation, the characteristic po- 
tentials it 7 (r) are determined by 



V 2 w Q (r) = 47re 2 



dr'n(r, r> a (r') 



dn(r, a) 
dE 



(59) 
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with the boundary conditions discussed in Section IIIII 
The above equation has the same structure as Eq. (fTB")) . 
(We postpone a detailed comparison between both equa- 
tions to the forthcoming Section.) Both terms in the 
r.h.s. of Eq. (f59|) stem from of the electronic charge im- 
balance n^(r)M- 

By using (|30|) and integrating by parts, the important 
relation 



E 



dn(r, a) 
dE 



= / dr'n(r,r') 
Jv 



(60) 



is obtained. This relation holds also beyond the flat and 
wide band approximation^ 

We now sum Eq. (f59|) over all leads a to write 



V 2 Vu Q (r) = 47re 2 / dr'n(r,r') 
J v 



.(61) 



Recalling the boundary conditions for u a we find that 
Y] a u a = 1, formally recovering one of the sum rules put 
forward by Biittiker 8 . This is a quite simple way to prove 
that, in the Hartree approximation, the NEGF formalism 
is manifestly gauge invariant. Conversely, the scattering 
approach uses J2 a Ma = 1 to obtain (fT3|) . 

Another way to picture the sum rule J2 a u a = 1 is by 
observing that it automatically guarantees that n^(r) 



remains invariant under the global shift V a — » V a + Vq 
(as shown to all orders in the previous section). 

As experimentally established 1 -^^ and theoretically 
discussed , 10 i 15 the Onsager-Casimir reciprocity relations 
do not hold for non-linear conductance. In the formal- 
ism we present, this is manifest in Eq. (|54[) . While 
5T a p/SU(r) computed at {V^,} = is even in magnetic 
field, in general u 7 (r, B) ^ u 7 (r, —B). The later can be 
seen from Eq. (|59|) : Albeit the Lindhard function is even 
in B 



(62) 



n(r,r';B) = IT(r, r'; —B), 



as a consequence of the "microreversibility" relation 1(501 
G r ^(r,r';E,B) = G r{a \r' , r; E, -B), in general 



dn(r,a,B) ^ dn(r,a,—B) 



dE 



dE 



(63) 



C. Arbitrary order 

Recent experiments measured higher order conduc- 
tance coefficients^, calling for a theoretical analysis of 
higher conductance coefficients. The general expression 
for G a p 1 ---i3 J in terms of the "generating functional" T a p 
is 



2e 2 [°° ( df \"^ {-iy+\ 

00 ^ ' 1=1 n=0 



^ 1 / * ■ ■ ■ d < ww) - ■ • m^$m suK) K ^ ' ' ' ■ r - } ' (64) 



where is defined as 



(65) 



{V T }=0 



and related to the characteristic potentials by (|3]). For n = 0, we define _K"(°) = <5j;. 

Let us explain the structure of Eq. (f64|) . We identify the term containing V to the power J in the I—V characteristics 
Eq. ([5]) with the conductance coefficient G a .p 1 ,---i3.j ■ The V J term comes from the product of the expansions of T a p 
and fo in (|28p. The first sum in Eq. run over i, the power of V that stem from the expansion of fo- The J — I 
derivatives of T a p with respect to V give raise to higher order characteristic potentials and functional derivatives of 
the kind 5T a p/8U{r'). The second sum in Eq. (|64[1 run over n, that represent the number of derivatives dy that 
become u~/(r')6/5U(r'). The remaining J — I — n derivatives are responsible for higher order characteristic potentials 
u^fo.... The expansion of fo is straightforward, but to factorize the (— 5b/o) term, we have to integrate by parts. As 
a result, we obtain additional I — 1 functional derivatives acting on T a p. 

Note that, as discussed in Sec. [II] the indices /3i • • • fii are not symmetrized. By combinatorial arguments it is not 
difficult to find the Taylor coefficients of the I — V expansion. 

Similar expressions have been obtained classically — The connection between the classical and quantum results is 
not clear yet. 

The expression for the characteristic potential of order J is obtained from the expansion of G < in Eq. 



k tcrr 



/7 jyi r\ r\ J 

o-7nT-~ 7nT E< r l G 5 v ^ G o ' ' ' v ^ G o E K E ~ eV ^ T <* 
2tt dv Pl dV 0J ^ ^ 
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J-k 



j terms 



x J2 G o ^ffGg • ■ • KffGg |r) 

j=o 



(66) 



{v.,}=o 



The effective potential V e g = eU — eU cq does not change upon a global shift V a — > V a + Vq. 

Equations ([64]) and ([66]) allow us to solve the non-linear problem to arbitrary order. We checked that, in lowest 
order, these equations lead to the results discussed in the previous sections, as they should. For the sake of illustration 
we explicitly show the third order conductance coefficient, namely, 



dE - 



dfo 



dE 
dri I dr 2 



ST af3 



f (ri)u s (r 2 ) - 5^us(ri) + -6^5 lS 



(67) 



and the corresponding equation for u a p 
V 2 u al3 (r) = -47re 2 



dn(r) d 2 n(r) d 2 n(r, 0) d 2 n(r,/3) , , d 2 n(r,a) . . 

V( r ) + e dE 2 u a (r)uf}{Y) + e—p^ — d a p - e — — u„(r) - e — — u^(r) 



dE 



dE 2 



dE 2 



dE 2 



(68) 

For simplicity, u a p is written in the Thomas-Fermi approximation. It is simple to explicitly show that ^ Q u a p = 0, 
respecting gauge invariance. 



r 



VI. CONNECTION WITH THE 
SCATTERING-MATRIX APPROACH 

We now establish the equivalence between the NEGF 
results for nonlinear elastic electronic transport and those 
obtained by the scattering approach, as formulated by 
Biittiker and collaborators^ and summarized in Section 

IE 

The underpinning elements of the scattering approach 
are: (a) The scattering matrix is viewed as a func- 
tion of the electron energy and a functional of the 
non-equilibrium electrostatic potential in the conductor. 
S(E, U(r)) is used to construct a generating functional to 
obtain the nonlinear conductance coefficients, (b) Phys- 
ical arguments are used to write a Poisson equation re- 
lating U(r) with source terms expressed as functions of 
S. Both S and U are solved self-consistently. 

We begin by writing the standard (equilibrium) reso- 
nance scattering matrix 4 ^ 

S ba (E) = S ba - 2iriY,(paPb) 1/2 V; b {G r Q (E)]^V„ a , (69) 



where a and b label propagating modes in the leads and 
the conductor Green's function is calculated for the equi- 
librium electrostatic potential U cq (r). As the bias is in- 
creased U(r) is modified, as well as the conductor reso- 
nances and the channels thresholds. Provided that the 
applied bias does not add new physical processes to the 
transport problem, like for instance phonons, the "equi- 
librium" S'-matrix can be generalized to 

S ba (E,U(r)) = S ba - 2 m Y J {PaPb) 1/2 V; b [G r {E)]^V va . 



(70) 



Here the non-equilibrium U (r) is contained in the con- 
ductor Green's function G r . By recalling the definition 
of the decay widths T , we obtain 44 



A a p = Tr 



= —T a p 



(71) 



and conclude that the scattering and the NEGF ap- 
proaches formally give identical expressions for the cur- 
rent. 

We now examine how U(r) is treated in both ap- 
proaches, namely, we compare Eqs. (fl~3|) and (|59|) . We 
begin by comparing the injectivities. By making explicit 
the sums over channels, the scattering approach injectiv- 
ity ([5]) becomes 



dn s (r, a) 
dE 



— E 

2m ^ 



00 dE 
.™ 2tt 



dfo 
dE 



E 

be/5 



ss, 



ba 



SS; 



eSU{r) eSU{r) 



Sba 



(72) 



evaluated at {V 7 } = 0. Equation (|70|) renders a quite 
amenable path to calculate the functional derivative 
SSl a /5U(r): First, we use G r <5[(G"') _1 ] + (SCT^GT)- 1 = 
to write 

5Sba 



eSU(r) 



= +2iri(p aPb ) 1 / 2 



E W(^)w 



8[G r {E)-\, v , 
eSU(r) 



[G r {E)] v , v V va . (73) 



Second, noting that (G r )~ 1 = E - H - eU + iT/2 and 
Uuv = J dr' ' {ii\y')U(y i ){y'\v) 1 we readily write 

5 



eSU(r) 



[G r (E)-% u = -{»\r){r\v). 



(74) 
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Plugging Eqs. (g9j and into (J72J) and using d74j) we 
show that dn s (r,a)/dE exactly coincides with Eq. ([55]). 
obtained within the wide and flat band approximation. 

Following the same steps, we also find that the scat- 
tering approach emissivity is given by 



dn(a, r) 
dE 



The Thomas-Fermi approximation is key to have a 
closed calculational scheme in terms of scattering matrix. 
As standard, it is assumed that u a (r) shows a slower co- 
ordinate dependence than II (r, r') to write 



f dr'n(r,r>a(r') w u a (r) / dr'II(r,r') 
Jv Jv 



u a {r) 



v 

dn(r) 
dE 



(76) 



where we use (|60|) . In this limit the NEGF Hartree equa- 
tion ([55)) reduces to the scattering Poisson equation ([T3")) . 
It is interesting to note that the diagonal part of II(r, r') 
can be written as 



n(r,r) 



-2i 



dE 
2^ 



dfo] 
dE ) 



(r\G a (E) - G r Q (E)\r) 



which, integrating by parts and using (|30p . gives 

dn(a, r) 



n(r,r) = ^- 



dE 



(77) 



VII. NEGF IMPLEMENTATION FOR THE 
HARTREE-FOCK APPROXIMATION. 

In order to improve our mean field description of non- 
linear conductance, we now include the exchange inter- 
action term, and consider the many-body problem in the 
Hartree- Fock approximation. The non-local nature of the 
exchange interaction prevents the description of the non- 
linear conductance in terms of (local) characteristic po- 
tentials defined in Eq. ([3]), making a description in terms 
of the scattering matrix approach unpractical. However, 
as before, it is possible to construct a self-consistent ex- 
pansion for the I — V characteristics, and treat the prob- 
lem by NEGF. 

Let us consider the interacting Hamiltonian 



where 



He = Ho + Hi 



(78) 



(79) 



describes the single-particle terms of conductor Hamil- 
tonian, and replaces the operator He addressed so far, 
and 

Wint = V ^ s 4«4 a ' < V d &' » ( 8 °) 



ss' fivyS 



with 



Note that Eq. |68J coincides with the corresponding one 
obtained in the scattering approach within the Thomas- 
Fermi approximation^ This suggests that, within the 
approximations discussed in this Section, both ap- 
proaches are equivalent to all orders. 



V Min s = J dridr 2 0*(ri)0*(r 2 )F(ri - r 2 )0 7 (r 2 )05(ri) . 

(81) 

In the Hartree-Fock approximation, the interaction 
term reads 23 



H, 



HF 



^W? ( d ls d Ss)di s ,d ls , - (dl s d ls ,)di s ,d Ss + ^(4* d ««)(4 a ' < V) - l( d U d ^')( d ts' d Ss) 



(82) 



Since the H^f is a bilinear operator, it is straightforward to obtain 



G r< - a \E) = [(G r W(£J))- x - H 



HF 
Int 



(83) 



where the H^f matrix elements are 



(84) 



The lesser Green's function is given by and the self consistent equations reads 

G<<*> - G r (E)T: < (E)G a (E) (85, _ _. j ^ ^ _ _^ j gg;^. 

(86) 
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Gauge invariance is shown to hold by the same arguments 
used in Section [TV] 

Eqs. (|83| . (|85p and (|86f provide the elements to write 
a power expansion of I a , in analogy to Section |Vj The 
non-local nature of the exchange interaction is encoded 
in Eqs. (f81"j) and (|82|) . and leads to a more involved self- 
consist scheme than that of Section HVl 

VIII. CONCLUSIONS 

We studied the nonlinear phase coherent quantum 
electronic transport properties of nanoscopic devices us- 
ing the Noncquilibrium Green's function method. This 
method allows us to express the I — V characteristics of 
a given system to arbitrary powers of the applied volt- 
ages in terms of equilibrium Green's functions. We show 
that the formalism is gauge invariant, provided that U (r) 
is calculated self-consistently and the induced charge is 
well localized. The latter condition is key to partition the 
system as in Eq. (fl4l) . the starting point of our discussion. 

We explicitly establish a connection between the 
NEGF method and the scattering approach. This is 
done by analyzing the first nonlinear contributions to 

(2) (2) 

the current, namely, I a ■ We show that I a obtained 
by NEGF at the Hartree level reduces to the scatter- 
ing matrix result in the Thomas- Fermi limit (and by us- 
ing the wide band approximation). It should be noted 
that while in the scattering approach gauge invariance is 
used to construct the Poisson equation, NEGF renders 



gauge invariance automatically. These observations sug- 
gest that NEGF provides a framework for treating the 
many-body problem at a more accurate level of approxi- 
mation. In Scc. lVIIl we discuss the Hartree-Fock approx- 
imation, very amenable to treat with NEGF, but clearly 
unsuited to the scattering approach, which is restricted 
to local potentials. 

We also analyze the electronic transport symmetry 
with respect to magnetic field inversion. In particular, we 
discuss the consequences of "microreversibility" for the 
conductance coefficients Gap... using the NEGF method 
for the second order coefficients. The general conclusions 
are the same as the ones obtained from the scattering 
approach. We then generalize to the theory arbitrary or- 
der in V, which is useful to address nonlinear transport 
experiments, such as Ref. 0. 

In general, as the bias is increased, very quickly in- 
elastic channels are opened. The inclusion of inelastic 
processes in the formalism and the development of ap- 
proximation schemes to solve such problem is one of the 
next main goals to pursue. 
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